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Abstract 


Spar's  [11]  vertically  integrated,  diabatic  model  is 
modified  so  that  it  possesses  an  invariant  mass  integral  of 
the  sum  of  enthalpy  and  the  kinetic  energy  of  the  non-divergent 
wind  in  the  adiabatic  case.  The  model's  linear  characteristics 
are  discussed,  and  a  series  of  numerical  experiments  utilizing 
idealized  initial  conditions  is  proposed. 


Introduction 


The  results  [12]  obtained  with  the  diabatic  model  derived  by 
Spar  [11]  indicated  that  diabatic  processes  were  capable  of  influencing 
the  development  of  a  cyclone.  The  predictions,  however,  were  found 
to  be  subject  to  a  spurious  anticyclogenesis.  This  effect  was  attributed 
to  an  inconsistent  treatment  of  the  two  components  of  the  curl  of  the 
vertical  advection  of  the  horizontal  wind  term  in  the  vorticity  equation. 
Arnason  and  Carstensen  [2]  had  previously  discussed  such  an  effect 
but  it  was  considered  to  be  important  only  in  hemispheric  scale 
predictions.  As  a  result  of  the  apparent  significance  of  inconsistencies 
of  this  type  for  even  short-range  predictions,  we  were  led  to  consider 
other  possible  inconsistencies  in  the  formulation  of  the  model.  Three 
types  of  inconsistency  were  found  to  exist.  The  first  involves 
assumptions  regarding  the  modeling  of  the  vertical  distribution  of 
various  model  parameters.  For  example,  the  assumed  vertical  tem¬ 
perature  profile  is  taken  as  invariant.  This  assumption  is  in  general 
not  consistent  with  the  predictive  equations  of  the  model.  It  was  con¬ 
cluded,  however,  that  this  kind  of  inconsistency  is  inherent  in  all  but 
the  simplest  atmospheric  models.  Even  a  model  with  many  infor¬ 
mation  levels  requires  the  implicit  assumption  of  a  prescribed  vertical 
distribution  within  the  layers  separating  the  information  levels,  and 
this  distribution  (e.  g.  linear  variation)  is  generally  incompatible 
with  the  predictive  equations.  It  was  decided  that  no  attempt  would 
be  made  at  this  time  to  eliminate  the  first  type  of  inconsistency  on 
the  assumption  that  it  does  not  lead  to  serious  errors. 
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Secondly,  following  the  ideas  of  Lorenz  [9  |,  we  found  that  the 
model  equations  possessed  an  inconsistent  energy  integral.  (Spar  [11] 
had  previously  discussed  the  inconsistency  in  the  circulation  integral.) 

It  was  concluded  that,  although  such  an  inconsistent  energy  integral 
might  not  be  crucial  for  an  adiabatic  model,  it  would  deprive  a  diabatic 
model  of  an  important  balance.  If  spurious  energy  sources  are  present 
in  an  adiabatic  version  of  a  model,  it  is  likely  that  the  diabatic  model 
will  also  contain  spurious  sinks  and  sources,  but  their  nature  cannot 
be  anticipated.  The  major  part  of  this  paper  is  concerned  with  this 
problem,  and  presents  a  simplified  version  of  the  diabatic  model 
proposed  by  Spar  [11]  which  has  been  modified  to  possess  a  consistent 
energy  integral.  The  model  is  derived  and  formulated  as  a  set  of 
difference  equations  which  are  to  be  integrated  numerically  using  a 
series  of  idealized  initial  conditions.  It  is  at  this  point  that  we  come 
to  the  third  type  of  inconsistency. 

Generally,  the  discussion  of  the  merits  of  difference 
schemes  has  been  based  on  the  order  of  magnitude  of  their  truncation 
errors.  A  different  but  related  problem  is  connected  with  the 
possibility  that,  although  the  differential  equations  of  the  model  may 
possess  a  consistent  energy  integral,  the  finite  difference  equations 
may  not  preserve  this  property.  Shuman  [10]  has  devoted  considerable 
effort  to  the  study  of  this  problem  utilizing  the  so-called  primitive 
equations.  He  found  two  schemes  for  approximating  the  non-linear 
advection  terms  which  preserved  the  energy  integral  of  his  fundamental 
equations.  However,  no  comparable  work  has  been  reported  regarding 
the  proper  choice  of  difference  operators  for  terms  appearing  in  the 
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filtered  equations.  Gates  [8]  has  recently  reported  a  related  difficulty 
encountered  in  the  hemispheric  integration  of  an  energetically  consistent 
two  layer  model  with  variable  static  stability.  In  an  effort  to  minimize 
this  error,  we  have  incorporated  in  the  difference  equations  an  approxi¬ 
mation  to  the  Jacobian  operator  suggested  by  A.  Arakawa  [1]. 


1.  Some  integral  properties  of  the  quasi -static  equations 

For  an  ideal,  inviscid  gas  the  hydrodynamic  equations  in  the 
p-coordinate  system  [5]  maybe  written 


36 

W 


^  +  V  •  q  W  +  Ik  •  V  X 


+  V-(q  Ik  X  \V)+ V 


(“ 


3\V 

3p 


+  V*V(^  +  K)  =  0 


-V 

9p 


W  s  -  6 


d<p  _  RT 
3p  ■  p 


(1.  1) 
(1.2) 

(1.3) , 

(1.4) 


(1.5) 


Equations  (1.  1)  and  (1.2  govern  the  relative  vorticity  (  C)  and  the 
divergence  (6)  of  the  horizontal  wind.  The  following  notations  are  used: 


q  =  t  +  f  ,  the  absolute  vorticity, 
f  ,  the  Coriolis  parameter, 

u  =  ^  ,  the  vertical  p -velocity, 

V  ,  the  horizontal  gradient  operator, 

and  Ik  ,  a  unit  vector  normal  to  the  x-y  plane  directed 

outward  from  the  earth. 


the  geopotential  pressure,  temperature,  gas 
constant,  specific  heat  at  constant  pressure,  and 
heat  supplied,  respectively, 

the  kinetic  energy  per  unit  mass  associated  with  W. 
Equation  (1.3)  is  the  equation  of  continuity,  and  equation  (1.5)  is  the 
first  law  of  thermodynamics  for  an  ideal  gas. 

We  propose  to  examine  two  integral  properties  of  the  set  (1.1) 
through  (1.5)  for  the  purpose  of  establishing  the  consistency  of  simplify¬ 
ing  approximations  which  may  be  introduced.  The  paper  by  Lorenz  [9] 
serves  as  the  guide  for  most  of  this  analysis. 

In  what  follows  we  shall  have  frequent  recourse  to  the  diverg¬ 
ence  theorem, 


A  B 


where  dcr  is  an  element  of  surface  area,  B  is  the  boundary  curve  of 
the  area  A,  and  df  is  line  element  of  B.  94^/ 9n  is  the  derivative 
of  4*  in  the  direction  of  the  outward  normal  to  B,  and  4  and  f  are 
arbitrary  scalars. 

The  following  symbols  will  be  employed. 

(~)  =  /  (  )  d(r  ( 

A 

_  Po 

()=/()  dp  ( 

Pt  ■ 

In  (1.8)  p^  and  p^  are  fixed  values  of  the  pressure,  p,  representing 
the  upper  and  lower  boundaries  of  the  atmosphere. 

QC 

If  the  operator,  (  )  ,  is  applied  to  equation  (1.1)  we  obtain. 


<p,  p,  T,  R,  Cp,  H, 


and  K, 
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m 


|t+^Ti-7(|k  X  V4')di  +^T,  ||di 
B  B 


B  B 

The  velocity  W  has  been  represented  above  by  the  sum  of.  a  solenoidal 
and  an  irrotational  vector.  Thus 

W  =  !k  X  V4^  +  Vx  >  (1.10) 


where  is  a  streamfunction  and  x  is  a-  velocity  potential.  If  the  hori¬ 
zontal  wind  has  no  component  normal  to  B,  and  if  w  vanishes  on  B,  then 


(1.11) 


Whenever  the  region  A  is  such  that  the  line  integrals  on  B  vanish,  we 
shall  refer  to  A  as  a  "closed"  region.  The  result  (1.11)  is  one  of  the 
integral  properties  of  the  system  of  quasi -static  equations. 

For  the  formulation  of  the  energy  integral  it  is  convenient  to 
partition  the  kinetic  energy  of  the  horizontal  wind  into  components, 

K  s  +K^  (1.13) 

in  which 


and 


^2  “  ^ 
K3  ®  Vx  • 


Hence 

|K,y»,v|+.Vx-v|ft^(Vx-lkx  V+1. 


(1.13) 

(1.14) 

(1.15) 

(1.16) 
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^  ^  ^  If] 

The  right  hand  side  of  (1.17)  may  be  evaluated  by  using  equations  (1. 
and  (1.2).  Under  the  assumption  that  the  line  integrals  vanish, 


qW.V4'  +  Ik 


(vt  X §2) 


From  (1.2), 


X  |^  =  nlkx  W'Vx  +  7x’VK  +  wVy 

From  the  continuity  equation. 


Vx-VK  =K 


8w 

8p 


It  follows  that 


8K 

5F 


From  the  thermodynamic  energy  equation  (1.5), 


8c  X  w  _>  ^ 

•  Vx  +  ^  (CpT  +  0)W 


Addition  of  (1.22)  to  (1.21)  leads  to 
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J 


5^  (K  +  CpT  +  0)w 


Thus,  in  the  adiabatic  case  the  total  energy  of  the  system  is  con¬ 
served  if  u  vanishes  on  the  upper  and  lower  boundaries. 

On  the  assumption  that  the  vertical  velocity,  w,  vanishes 
at  p^  and  p^  ,  (1.21)  and  (1.22)  yield 


(1.17) 

1) 

(1.18) 

(1.19) 

(1.20) 

(1.21) 

(1.22) 

(1.23) 
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and 


-  VX'V0 


8c  T 


(1.24) 

(1.25) 


Vx  •  W  = 


(1.26) 


In  view  of  (1.24)  and  (1.25),  (1.26)  implies  that  a  positive  conversion 

of  potential  energy  into  kinetic  energy  requires  that  a  negative  cor¬ 
relation  exist  between  oj  and  the  specific  volume.  Such  a  correlation 
exists  when  warm  air  ascends  and  cold  air  sinks  on  the  average. 

The  two  integral  theorems,  derived  above  as  equations  (1.11) 
and  equations  (1.24)  and  (1.25)  are  regarded  as  fundamental  properties 
of  the  quasi-static  system.  In  what  follows  certain  simplifying  ap¬ 
proximations  will  be  introduced  into  the  quasi -static  system.  We  shall 
judge  the  validity  of  these  approximations  by  seeing  how  well  they  pre¬ 
serve  the  integral  theorems. 

It  is  seen  that  if  the  kinetic  energy  of  the  irrotational  part 
of  the  wind  is  deleted  from  the  total  kinetic  energy  in  eq.  (1.24),  the 
integral  theorems  will  hold  for  the  general  balanced  quasi -static 
system. 

Thompson  (14]  demonstrated  the  necessity  and  sufficiency 
of  this  approximation  as  a  means  of  eliminating  internal  gravity 
oscillations  in  the  system.  However,  the  general  balanced,  quasi¬ 
static  system  is  too  complex  to  be  attractive  as  a  basis  for  an 
atmospheric  model.  In  order  to  discuss  additional  simplifications, 
it  is  convenient  to  expand  the  modified  version  of  (1.2), 
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-  V'  riV'i'  +  V  •  -  V  •  VKj 
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8 


V  •  |k  X  Vy  +  V’  u  X  V4^  +  V*  w  Vy  +  V*  ^1^2  +  V'VKji  -  0 

(1)  (2)  ^(3)  (4)  {5)-* 


(1.27) 


The  bracketed  terms  above  have  usually  been  neglected,  while  the 
remaining  three  terms  are  considered  to  comprise  the  balanced  wind 
equation. 

We  wish  to  determine  how  the  neglect  ot  the  bracketed  terms 
affects  the  integral  theorems.  For  this  purpose  we  have  numbered  the 
various  terms. 

Consider  first  the  energy  integral.  If  we  premultiply  each  of 
the  bracketed  terms  in  (1.27)  by  the  velocity  potential  and  apply  the 
operator,  (  )  ,  we  obtain. 


0  -  wVx  •  0^  |k  X  V4'  •  w  -  Vy  '  VK2  ■  '^y  •  VK^ 

11)  (2)  (3)  (4)  (5)" 


(1.28) 


The  terms  numbered  (3)  and  (4)  in  equation  (1.28),  may  be  combined 
by  use  of  the  continuity  equation.  If- we  apply  the  operator,  (  ),  to  the 
result  we  find  that  this  pair  vanishes.  Consequently  the  neglect  of  the 
terms  numbered  (1),  (3)  and  (4)  in  equation  (1.27),  does  not  require 
any  modification  of  the  vorticity  equation  in  order  to  maintain  the 
energy  integral  property  of  the  general  balanced,  quasi -static  system. 
It  can  be  shown  that  the  neglect  of  the  terms  numbered  (2)  and  (5)  does 
require  that  the  vertical  derivative  appearing  in  the  vorticity  equation 
be  modified  by  deleting  the  irrotational  wind  component. 

Thus,  if  we  rewrite  the  vorticity  equation  in  the  form. 
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8C 

5r 


+  V-  tii’k  xVi)-'  +  V*  ■nVx  +  Ik*  V  X 


k  X  ViJj 


=  0 


the  balanced  wind  equation  in  the  form 


(1.29) 


V*  +  V*V0  +  7*VK^  =  0 


will  be  consistent  with  the  energy  integral 


=  Vx*W 


(1.30) 


(1.31) 


The  modified  vorticity  equation  (1.29),  still  possesses  the  integral 
property,  (1.  11). 

A  further  simplification  which  has  been  widely  used  arises 
from  the  desire  to  eliminate  the  complexities  arising  from  the  last 
term  in  (1.29).  Arnason  and  Carstensen  [2]  have  shown  that  incon¬ 
sistencies  in  the  numerical  approximation  of  this  term  can  lead  to 
deleterious  results  in  forecasting  over  a  large  area.  It  is  readily 
seen  that  if  this  term  is  neglected  in  (1.29)  then  for  energetic  con¬ 
sistency  the  term,  V-VKj^  ,  should  be  eliminated  in  (1.30).  It  is 

probably  easier  to  justify  this  modification  of  the  system  by  appealing 

2  2 

to  the  smallness  of  V  in  comparison  with  V  <p  within  the  context 
of  equation  (1.30).  Nonetheless,  the  system. 


+  V*  (^Ik  X V4j)  +  V*  (tiVx)  =  0 

V.(tiV+)  = 


(1.32) 


constitutes  an  energetically  closed  system  with  the  same  integral 
properties  as  the  general  balanced,  quasi-static  system. 

Because  the  balance  equation  of  system  (1.32)  is  non-linear 
in  4’  and  because  the  relative  vorticity  is  frequently  small  compared 


with  the  absolute  vorticity,  many  numerical  prediction  models  have 
been  based  on  the  quasi -geostrophic,  balance  equation, 

(1.33) 

c 

in  which  is  given  a  constant,  mean  value.  The  form  of  the  vorticity 
equation  which  is  energetically  consistent  with  (1.32)  is, 

+ V-(iiik  xV4^)  +  f^V^X  =  0  (1.34) 

This  set  also  preserves  the  areal  mean  value  of  the  relative  vorticity. 

One  further  simplification  is  useful  in  this  quasi -geostrophic 
system.  It  is  not  necessary  for  consistency  with  the  energy  or  vorti¬ 
city  integrals,  but  it  is  probably  consistent  with  the  tangent  plane 
approximation.  This  simplification  involves  writing  the  second  term 
in  (1.34)  as 

V-(Tilk  xVi)^)  -  J(4',4)+  (1-35) 

in  which  J(4',  C)  is  the  Jacobian  operator  and  p^  is  a  constant  mean 
value  for  the  gradient  of  the  Coriolis  parameter. 

In  the  preceding  analyses  we  have  not  found  it  necessary  to 
introduce  the  first  law  of  thermodynamics.  We  have  implied,  however, 
that  the  horizontal  wind  appearing  in  the  equation  was  modified  as  a 
result  of  our  modification  of  the  divergence  equation. 

The  set  of  equations  that  has  been  chosen  as  the  basis  for 
the  prediction  model  is  the  following: 
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If  ♦  *  Pr  K  +  ”  » 


8oj  ^2 

§?=  X 


and 


^  =  -  jBlZ 

8p  ‘  p 


ax 

■§t 


+  V- 


\v+ 


=  v-^-\v+  ii 

p/  p  p 


(1.36) 

(1.57) 

(1.38) 

(1-39) 

(1.40) 


Although  this  set  is  considerably  simpler  than  the  general  quasi-static 
equations,  the  fact  that  it  possesses  the  integral  properties  discussed 
above  recommends  it  as  a  useful  system  for  studying  certain  simpler 
aspects  of  the  transformation  of  energy  and  growth  of  circulation 
within  the  atmosphere. 


2.  Derivation  of  the  integrated  model 

In  the  preceding  section  we  considered  several  simplified 
versions  of  the  hydrodynamic  equations.  A  simple  system  that 
possesses  an  appropriate  energy  integral  was  adopted  for  use  in  this 
study.  We  shall  now  discuss  a  model  formulation  of  these  equations 
in  which  we  follow  in  many  essential  respects  the  approach  used  by 
Spar  [11].  This  is  a  particular  type  of  the  thermotropic -model  de - 
veioped  by  Thompson  [13].  (Another  such  model,  which  is  somewhat 
more  general  than  Spar's,  was  used  by  Berkofsky  [4],  but  no  numeri¬ 
cal  experiments  with  that  model  have  been  reported.  ) 


The  thermotropir  models  are  based  on  the  assumption  that  the 
isobaric  temperature  gradient  does  not  change  in  direction  but  may  var 
in  magnitude  as  pressure  changes.  One  may  express  this  symbolically 
as 

V  T(x,  y,  p,  t)  =  A(p)V  T(x,  y,  t)  (2.1) 

in  which  we  have  used  the  "bar -operator",  defined  by 

.  Pt 

(  f  (  (2.2) 

Po  -  Pt  p^ 

Since  this  operator  will  find  frequent  use  later  on,  it  is  advantageous 
to  digress  for  a  moment  to  discuss  it.  If  we  imagine  the  portion  of  the 
atmosphere  for  which  we  wish  to  carry  out  our  predictions  to  be  con¬ 
tained  between  two  pressure  surfaces  p^  and  p^  (p^  ^  Pj,)>  then  the 
mean  value  of  any  property  of  the  model  atmosphere  may  be  defined  by 
(2.2).  In  this  model  we  shall  use  p^  =  1000  mb  and  p^  =  200  mb,  and 
also  assume  that  these  pressure  boundaries  are  rigid.  The  latter 
assumption  implies  that  co  (=  dp/dt)  vanishes  at  p^  and  p^  .  The  value 
of  any  parameter  at  p  =;  p^  will  be  designated  by  a  subscript,  o. 

The  difference  between  the  mean  value  of  a  parameter  and  its  value  at 
p  =  p^  will  be  denoted  by  a  subscript,  T. 

Returning  to  the  thermotropic  type  model,  we  may  observe  that 
the  variation  of  temperature  prescribed  by  (2.1)  is  sufficient  to  allow 
a  simple  form  of  baroclinic  development.  By  virtue  of  (2.1)  we  have 
upon  integration,  the  general  temperature  distribution, 

T(x,  y,  p,  t)  =  A{p)  T(x,  y,  t)  +  F(x,  y,  t)  .  (2.  J) 

Following  Spar  we  set  F(x,  y,  t)  =  0,  The  specification  of  A(p)  m.ay 


be  made  in  the  form 
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A(p)  =  a(p/p^)’"  (2.4) 

in  which  a  and  b  are  arbitrary,  but  not  independent.  Using  the  NACA 
standard  atmosphere  as  a  guide,  the  parameter,  b,  was  given  the 
value,  0.200.  The  requirement  that  the  mean  value  of  A(p)  be  unity 
may  then  be  used  to  specify  a,  which  was  found  to  be  1,123. 

The  vorticity  equation  is 

(2.5) 


in  which  is  the  relative  vorticity  of  the  horizontal  velocity,  4*  is 

the  streamfunction  of  the  horizontal  velocity,  (3^  (1.45  X  10  cm  ^sec 

is  a  constant  value  assigned  to  the  variation  of  the  Coriolis  parameter, 
-4  -1 

f^  (10  sec  )  is  a  constant  value  assigned  to  the  Coriolis  parameter, 
and  X  is  the  velocity  potential  of  the  horizontal  velocity. 

The  divergence  equation  is  replaced  by  an  equation  of  balance 
which  is 

f  V^4  =  (2.6) 

c 

in  which  <p  {  =  g^z  )  is  the  geopotential.  In  what  follows  we  will  use 
(2.6)  in  integrated  form,  viz.  , 

f^4  =  ?.  (2.7) 

The  hydrostatic  equation  for  an  ideal  gas  iTiay  be  written  as 


90 

9p 


RT 

P 


(2.8) 


in  which  T  is  the  temperature,  p  the  pressure,  a  the  specific  volume, 
and  R  is  the  gas  constant  for  dry  air  for  which  we  will  use.  the  value 
2.87  X  10^  ergs  gm  ^  deg  *. 


The  continuity  equation  is 


du 

in  which  w  is  the  vertical  p-velocity  (  ) 

The  first  law  of  thermodynamics  is 


^  =-J(4i.T)  -V 


•(t+  ^Vx  -  /-fT+  Vx-V  A+  ^ 


c  c 
P  P 


(2.9) 


(2.10) 


in  which  is  the  specific  heat  at  constant  pressure  (c^  -  1.00  x  10 
ergs  gm  ^  deg  ^),  and  H  is  the  rate  at  which  the  enthalpy  is  being 
changed  by  diabatic  processes. 

We  shall  also  require  equations  to  specify  the  value  of  H  but 
these  are  best  postponed  until  after  the  integrated  model  equations 
have  been  derived. 

By  using  (2.4)  and  (2.3)  in  the  hydrostatic  equation,  (2.8),  we 
may  arrive  at  the  model  representation  of  the  geopotential,  viz,  , 

0(x,  y,  p,  t)=  ^  (x,  y,  t)  +  E(p)  ^.p(x,  y,  t)  ,  (2.11) 

in  which 


E(p)  = 


1  -  A(p) 
a  -  1 


(2.12) 


From  (2.7)  and  (2.12)  ,  it  follows  that 


=  <t>  . 


^c^T  “  "^T  ’ 


and 


(2.11) 

(2.14) 


4'(x,  y,  t,  t)  =  i|>(x,  y,  t)  +  E(p)  ^^^.(x,  y,  t)  . 


(2.15) 


From  (2.14),  (2.11),  (2.12),  and  (2.8),  we  may  derive  the 
relationship. 
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+^=mT  (2.16) 

in  which 

'  _(a-t)R  2  -1,  -1  I  ->  I  ■7\ 

m  =  i — gyi —  =  1.76  X  10  cm  sec  deg  .  (2.17) 

c 

The  results  outlined  above  are  as  far  as  we  can  go  in  modeling 
the  parameters  of  the  equations  (2.5)  through  (2.10).  We  therefore 
make  an  additional  assumption  which  here  proves  to  be  convenient, 
and  in  our  analysis  of  the  model's  energetics  proves  to  be  quite 
essential.  It  consists  of  the  assumption  that  the  velocity  potential,  x> 
possesses  a  distribution  analogous  to  that  of  the  streamfunction,  viz.  , 

X(x,  y,  t,  t)  =  x{x,  y,t)  +  E(p)Xj(x,  y,  t)  .  (2.18) 

With  (2.18),  the  continuity  equation,  (2.9),  and  the  boundary  conditions 
we  may  derive  the  relations. 


X  =  0 

(2.19) 

w  =  F(p)  w  , 

(2.20) 

V  Xq.  =  -  V  W 

(2.21) 

v.{ 

f  E(p)  dp] 
^Po 

jj-  =  6,33  X  10’^  mb'^  , 

(2.22) 

F(p)  = 

P 

V  I  E(p)  dp  . 

(2.23) 

Po 

The  mean  vorticity  equation  is  derived  by  applying  the  "bar 
operator"  to  equation  (2.5).  The  result  is 

^  f  =  -  J{^,  f)  -  ^  ■  Pc  sl  ’ 


(2.24) 


in  which 


(2.25) 


E"  =  0.480 

It  should  be  noted  that  this  equation  for  the  mean  vorticity  differs  from 
a  barotropic  vorticity  equation  only  in  the  appearance  of  the  term 
involving  the  advection  of  shear  vorticity  by  the  shear  wind.  It  is  this 
term  alone  that  can  lead  to  development  of  the  mean  circulation. 

The  next  model  equation  may  be  derived  by  applying  the  vorticity 
equation  at  p  =  p^,  expressing  the  value  of  the  parameters  at  p^  in 
terms  of  the  mean  and  shear  parameters  and  then  subtracting  the  result 
from  the  mean  vorticity  equation.  One  obtains  in  this  manner  the  "shear" 
vorticity  equation, 

Birp  —2  _  _  S+rp  2 

-gf  =  (1  -  E^)  J(vp,p.;,p)  -J(4'.5t)  -Pc  -glT'^c^ 

By  using  equation  (2.21),  the  model  continuity  equation,  in  (2.26) 
we  may  replace  it  by 

-gi  =(1 -E^)  J{^,j.,t,p)-J(vP.t^)-J('4>,p,t)  -(3^  -g^  +f^Yw.  (2.27) 

We  may  now  consider  the  first  law  of  thermodynamics, 
equation  (2.  10).  Upon  application  of  the  bar -operator  to  (2.10),  the 
equation  reduces  to  the  form, 

^  =  -  J(^,f)  -  EA  J{v|i,p,  f )  -  r{p)  E(p)  V^Xt 

+  E^7Xp,*V^+ .  (2.28) 

In  (2.28)  we  have  introduced  an  approximation,  viz.  ,  that 


the  term 
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(t  +  r(p) 

\  py 

is  a  function  of  pressure  alone.  This  assumption  implies  that  the  static 
stability  is  a  function  of  pressure  alone.  In  view  of  (2.16),  the  second 
Jacobian  in  (2.28)  vanishes.  We  may  formally  introduce  (2.16),  (2.14), 
and  (2.21)  into  (2.28)  to  write 

_  _  ^f  m 

=  -  J(4^>.j,)  +  Er  Ymw+  —  H,  (2.29) 

P  P 

in  which  the  parameters  have  the  following  numerical  values: 

ET  =  8.69  'K 
Y  =:  6.33  X  10'^  mb‘^ 

,  ,.10  2  -1  o  -1 

m  =  1.76  X  10  cm  sec  deg 

f  =  lO"^  sec'^ 
c 

?  =  0.480 

7  2-2  -1 

c  =  1.00  X  10  cm  sec  deg 
P 

We  note  that  both  equations  (2.29)  and  (2.27)  involve  the  tendency 
of  the  shear  streamfunction.  By  taking  the  Laplacian  of  (2.29)  we  may 
combine  the  result  with  (2.27)  to  derive  an  alternative  diagnostic 
equation.  We  regard  the  resulting  equation  as  governing  the  value  of 
w  required  to  yield  a  flow  which  preserves  both  the  balanced  stream- 
function  and  the  diabatic  temperature  fields.  The  oo  equation  may 
be  written  after  some  rearrangement  in  the  form. 
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-J— rv^(J(4i.4<  ))  +  (!  -?) 

\Erm>  ErvmL  ^  1  T  1 


J(4^T.t)  -P. 


^T’*'  TT 


^'c  ^ 

-  V^(V+^-VXt,)  -  — ^ - 

EFvc^  EFvc 

P  ^  P 


(2.30) 


Equations  (2.30),  (2.26),  (2.24),  and  (2.21)  are  the  basic  equations 
of  the  integrated  model.  The  dependent  variables  are  4^,  4^^,  w,  and  X-j  • 
It  is  necessary  to  close  the  system  by  giving  auxiliary  equations  for  the 
specification  of  H  in  the  w-equation.  It  should  be  noted  that  only  two 
of  the  four  equations  are  prognostic;  the  other  two  equations  are 
diagnostic.  Given  the  initial  and  boundary  values  for  4'  and  4^.^.  , 
w  and  X-p  may  be  computed  from  the  diagnostic  equations  by  simul¬ 
taneous  iteration,  subject  to  prescribed  boundary  conditions.  One  may 
then  use  the  solutions  for  w  and  Xj  to  extrapolate  the  streamfunctions 
into  the  future  by  means  of  the  prognostic  equations.  Before  proceeding 
to  a  discussion  of  the  method  employed  for  computing  H,  we  shall 
examine  the  integral  properties  of  the  equations  derived  above. 


3.  Integral  properties  of  the  model  equations 

It  was  shown  in  section  1  that  the  differential  equations  upon 
which  the  integrated  model  is  based  possess  an  energy  integral,  com¬ 
patible  with  that  for  the  general  quasi -static  equations.  We  shall  now 
demonstrate  that  the  integrated  model  equations  also  possess  a  suitable 
energy  invariant.  Since  the  parameters  in  the  equations  of  the  integrated 
model  do  not  vary  in  the  vertical  we  shall  need  only  an  integral  operator, 
{  ),  which  we  define  as  a  surface  integral  in  the  x-y  plane  over  a 
"closed?'  region. 


The  kinetic  energy  of  the  non-divergent  wind,  K,  may  be 


partitioned  as  follows 


K  =  K,j,  +  E(p)K 

in  which 

- 2 - 

V4> 

^ - 2 - ^ 

and  K  =  V4’’ 


(3.  ]} 


(3.2) 

(3.3) 

(3.4) 


It  follows  from  the  nature  of  the  modeling  function,  E(p),  that  the  inte¬ 
grated  kinetic  energy,  K,  is  given  by 

Upon  introducing  (3,2)  and  (3.3)  into  (3.5),  the  local  time  derivative  of 
K  may  be  expressed  as 


f  =Vf-v|f+E 

When  the  areal  average  of  (3.6)  is  taken  and  the 
employed,  we  find  that 


9K 

at 


divergence 


(3.6) 

theorem 


(3.7) 


The  right  hand  side  of  (3.  7)  may  be  evaluated  by  use  of  the  vorticity 
equations,  (2.24)  and  (2.27).  We  obtain 


2  0 


and 


+T 


Using  (3.  8)  and  (3.  9)  in  (3.  7)  we  have  the  result, 


9R  _  ^ 


(3.  8) 


(3.9) 


(3.10) 


The  preceding  result  demonstrates  that  the  mean  value  of  the 
kinetic  energy  of  the  non -divergent  wind  is  not  conserved  within  the 
model.  Since  the  production  term  arises  in  the  shear  vorticity  equation, 
we  may  speculate  that  the  increase  in  kinetic  energy  appears  in  the 
shear  wind  first  and  is  redistributed  to  the  mean  wind  via  the  develop¬ 
ment  term  in  the  mean  vorticity  equation. 

We  may  now  consider  the  first  law  of  thermodynamics  governing 
the  integrated  value  of  the  enthalpy,  c  T.  From  (2.16)  we  may  write 


Thus , 


c  T  = 
P 


T 


iLc 

9t  p  m9t^T 


(3.11) 


(3.  12) 


The  right  side  of  (3.  12)  is  readily  derived  from  equation  (2.29).  We 
obtain  the  result, 

~'2 

It  should  be  noted  in  passing  that  the  appearance  of  the  factor  E 


in  the  first  term  on  the  right  of  (3.  13)  is  dependent  upon  our  having 
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chosen  the  same  modeling  function,  £l(p))  for  both  the  non -divergent 
and  irrotational  parts  of  the  horizontal  wind.  Upon  combining  (3. 10) 
and  (3.  13),  we  find  that 

a(K  +  C  f )  _ 

which  may  be  compared  to  eqs.  (1.24)  and  (1.25). 

The  conservation  of  vorticity  may  also  be  shown  to  be  a  property 
of  the  equations  of  the  integrated  model.  The  mean  square  value  of  the 
vorticity. 


(3.15) 


may  be  examined  too.  It  follows  from  (3.15)  and  the  vorticity  equations 
that 

|?=2?f^Y^T  (3.16) 

Since  it  can  be  readily  shown  that 

?  Y  =  -  FE"  (3.  17) 

the  result  (3.  16)  compares  reasonably  with  that  for  the  simple  hydrostatic 
vorticity  equation  (2.  5),  viz.  , 

%-=  -2  f  ,  (3.  18) 

^  c  8p 

within  the  framework  of  the  integrated  model.  Since  by  (2.  16)  we 
may  write 

(3.19) 

Then  (3.  16)  may  be  interpreted  as  implying  that  the  circulation  is 


intensified  when  w  is  correlated  with  the  mean  temperature's  field  of 

2 

"lumpiness".  For  example,  if  cold  pockets  (V  T  >  0)  were  on  the 
average  associated  with  ascending  air  (w  <  0)  or  warm  pockets  with 
descending  air  then  the  circulation  would  be  weakened.  Thus  in  some 
sense  we  may  regard  the  occlusion  process  to  be  contained  within. the 
model. 


In  connection  with  the  approximation  introduced  in  equati.on 
(2.  28),  the  following  analysis  is  offered  as  justification  for  the  use  of 
this  simplification.  The  first  law  of  thermodynamics  for  an  ideal  gas 
in  quasi -static  equilibrium  maybe  written  in  the  p -coordinate  system, 
as 


in  which 


^  =  -  J(4'.  T)  -  Vx  •  VT  +  wS  +  ^  , 


>  0 


(3.  20) 


(3.21) 


is  related  to  the  static -stability  (a8/9p).  In  most  quasi -geostrophic 
models,  two  approximations  are  made;  viz.  , 

Vx-VT  =  0  .  (3.  22) 

S  =  S(p)  or  constant.  (3.  23) 

When  the  equation  incorporating  (3.  22)  and  (3.  23)  into  (3.  20)  is  inte¬ 
grated  over  a  closed  mass  of  the  atmosphere  one  finds  that  there  can 
be  no  change  in  the  average  enthalpy  for  adiabatic  flows. 

From  the  general  form  of  (3.  20)  the  proper  change  in 
enthalpy  for  adiabatic  flows  is  seen  to  be  given  by 
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^=-yX-VCpT+  wSCp  .  (3.24) 


Now  one  of  the  weaknesses  of  an  integrated  model,  or  other 
models  of  limited  vertical  resolution,  is  an  inability  to  specify  the 
term  S  with  good  accuracy.  If  we  consider  the  consequences  of  making 
assumption  (3.23),  but  retaining  the  term  Vy  •  VT  ,  we  may  derive 
by  use  of  the  continuity  equation,  the  result. 


8c  T 


(3.25) 


The  result  (3.  25)  will  reduce  to 


8c  T 


(3.26) 


which  is  the  appropriate  form  (see  (1.26)  ]  only  if  the  stratification 
is  neutral.  In  particular  for  the  model  temperature  profile  used  irt 
deriving  the  model  parameters  in  section  2,  (3.25)  leads  to 


=  0.71  cua 


(3.27) 


which  represents  an  underestimate  of  the  conversion  of  enthalpy. 

The  simplification  which  we  have  used  in  the  equation  (2.  28), 
while  analogous  to  the  procedure  discussed  above,  leads  to  the  proper 
estimate  of  enthalpy  conversion.  To  elaborate  on  the  procedure,  we 
may  first  transform  (3,20)  into  the  equivalent  form. 


+  Vx-VT  +  Vx-V^+ ^  . 

p  p 


(3.28) 


We  now  define 


and  observe  that 


Now  for  adiabatic  flow  we  obtain  from  (3.  28)  the  result 


ac  T 
P 

“St 


Vx  *  V0  =  a  to 


{  3.  29) 


(  3.30) 


(3.31) 


The  result  (3.  31)  makes  it  clear  that  the  conversion  of  enthalpy  within 
a  quasi -static  atmosphere  depends  upon  the  work  against  the  "pressure 
gradient"  done  by  the  ir rotational  wind. 

The  form  of  (3.  28)  suggests  that  this  process  may  be  evaluated 
within  the  framework  of  a  semi -geostrophic  model  with  weak  vertical 
resolution  by  treating  F  as  a  function  of  pressure  alone.  It  is  this 
procedure  that  we  adopted  by  our  assumption  in  (2.  28). 

The  choice  for  the  numerical  value  of  the  parameter,  , 
which  appears  in  (2,  28)  was  derived  from  the  mean  values  of  the 
pararheter  S  (=  -  9r/8p)  published  by  Gates  [7]  and  from  the  model 
function  for  the  vertical  velocity,  F(p).  To  do  this  we  made  use  of 
the  relations 

-  ET  =  r  ^  =  FF  w  ,  (3.32) 

fF"  =  (FF)'  -  F  =  FS  ,  (3.33 

and  the  model  continuity  equation. 
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which  led  to 

Er  =  +  i  Fs  >  0  .  (3.  35) 

The  average  winter  and  summer  values  of  FS  were  computed 

to  be 

5.53  X  10'^  'K  mb'^ 

and  4.  45  X  lO'^  'K  mb'^  , 

respectively.  We  chose  to  use  a  value  biased  toward  winter  values, 

—  -2-1 

viz.  ,  FS  =  5.50  X  10  *K  mb  .  Using  this  value  in  (3.  35)  we  obtained 

-3  -1 

the  model  value  (  y  =  6.33  X  10  mb  )  , 

Fe  =  8.69  ‘K  . 

4.  The  diabatic  term 

Since  the  differential  equations  governing  the  integrated  model 
possess  an  energy  invariant  of  a  general  type,  it  is  reasonable  to  em¬ 
ploy  the  model  in  a  study  of  the  energy  transformations  associated 
with  diabatic  processes.  In  particular,  we  have  formulated  approxi¬ 
mations  for  the  release  of  latent  heat  by  condensation  and  the  transfer 
of  sensible  heat  from  the  sea  to  the  air.  Because  of  the  simple  vertical 
structure  of  the  model,  we  will  be  forced  to  confine  ourselves  to  the 
influence  of  horizontal  variation  in  heating  upon  the  development  of 
circulation. 

We  have  made  use  of  the  relationship  given  by  Spar  [11]  for 
the  estimation  of  the  magnitude  of  the  sensible  heat  flux  at  the  air- 
sea  interface.  In  modified  units  we  may  write  Spar's  formula  as 

Q  =  6.01  X  lO"'^  V  (T  -  T  )  , 

s  o'  s  o 


(4.  1) 
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in  which  is  the  rate  of  temperature  change  in  degrees  per  second, 

is  the  surface  wind  speed  in  cm  sec  \  surface 

temperature,  and  is  the  temperature  of  the  surface  air.  To 

account  for  the  reduction  in  wind  speed  in  passing  from  1000  mb  to  the 

surface  we  replace  in  (4.  1)  by  80%  of  the  wind  speed  at  1000  mb, 

0.8  IW  I  .  The  formiila  then  becomes, 
o 

=  4.81  X  10'^  IW  I  (T  -  T  )  .  (4.  2) 

s  o  s  o 

Since  this  formula  was  derived  for  cases  in  which  the  heat  flux  was 

directed  from  sea  to  air,  we  arbitrarily  set  Q  =0  when  T  <  T 

s  so 

No  heat  flux  is  computed  over  the  land. 

In  order  to  establish  a  formula  for  the  rate  of  heating  due  to 
condensation,  we  must  consider  an  equation  for  the  conveyance  of 
water  vapor.  The  parameter  used  to  represent  the  water  vapor  is  the 
specific  humidity,  q.  A  continuity  equation  for  q  may  be  written 
neglecting  evaporation  as  follows, 

|a=  -  W.ql-VqVx  ,  (4.3) 

in  which  (8q/8t)^  is  the  local  rate  of  decrease  in  specific  humidity 
due  to  condensation.  Application  of  the  bar  operator,  (  )  ,  to  (4.  3) 
yields , 

=  ■  v-,vx  (*-4) 

If  we  assume  that  the  distribution  of  specific  humidity  may  be  written 
as, 

q(x,y,  t,t)  =  L{p)q{x,  y,t)  ,  (4.5) 


then  the  equation  becomes. 
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(4.6) 
c 

In  qualitative  agreement  with  the  results  of  Benton  and 
Estoque  [3],  a  reasonable  average  value  for  the  parameter,  LE  , 
has  been  found  to  be  -0.6.  Because  the  integrated  model  is  really 
incapable  of  accurately  describing  the  vertical  distribution  of  specific 
humidity,  it  was  decided  to  simplify  the  term,  V-q  Vy  •  by  replacing 
q  with  an  average  value,  q^  .  By  introducing  the  model  continuity 
equation  (2.  21),  and  the  approximation  just  discussed,  we  may  re¬ 
write  (4.  6)  as, 

=  -  J(^.q)  -  LE  J(v|j,j.,q)  +  LE  .  (4.7) 

Recalling  that  LE  =  -0,6  ,  we  may  note  that  the  term  involving  w 
will  act  as  a  source  of  q  when  w  <  0  (ascending  air)  and  as  a  sink 
of  q  when  w  >  0  (descending  air).  Now  if  the  air  column  is  locally 
saturated,  one  may  expect  that  the  Jacobian  terms  will  rarely  tend 
to  increase  q  above  its  local  saturation  value.  However,  if  w  is 
negative,  the  term  involving  w  'will  act  to  increase  q  above  its 
saturation  value.  We  are  therefore  led  to  assume  that. 


If  =  -  J(4>.  q) 


LE  J(vtJ,j.,^  -  LE  V-qVx^ 


LE  q^  yw  , 


(4.  8) 


if  CO  <  0  and  the  air  is  saturated. 

To  formulate  this  hypothesis  quasi -analytically  we  introduce 


the  parameter,  6,  defined  as 
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c» 

II 

o 

if 

w 

>  0 

or  q  <  qg 

5  =  1 

if 

(0 

o 

V 

and  q  >  Tg 

Using  (4.  8)  and  (4.  9)  in  (4.  7),  we  may  write 


.  J(iJj,q)  -  LE  +(1  -  6)LE  yw 

Now  if  6  =  1,  the  rate  of  change  in  mean  temperature  in  the 
column  may  be  obtained  by  multiplying  (4.  8)  by  the  ratio  of  the  latent 
heat  of  vaporization,  Lv,  to  the  specific  heat  at  constant  pressure, 
c  .  We  obtain 

P  _ 

Lv  LE  q  Y 

Q  =  - - — w 

'  c  c 

P 

Using  the  values. 


Lv  =  600  cal  gm 

=  0.239  cal  gm  ^  deg 


LE  =  -  0.6 

q  =  2,5  X  10“^ 
o 

Y  =  6,33  X  10  ^  mb 


we  may  write 

=  IT  [wj  , 

in  which  tt  =  23.83  x  10'^  deg  mb'^  . 

The  precipitation  rate,  P,  may  also  be  computed  from  (4.  8) 
through  the  .relation 


(P^  -  P,) 

gp, 


(4.9) 


(4.10) 


(4.11) 


(4. 12) 


(4.13) 


in  which  g  is  the  acceleration  of  gravity,  and  is  the  density  of 
liquid  water.  Using  (4.  8)  we  obtain 
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P  =  r  luil  (4.  14) 

-3  -1 

in  which  r  =  7.60  x  10  cm  mb 

In  order  to  compute  the  saturation  specific  humidity  from  the 
model  parameters,  we  adopt  the  approximation  by  Spar, 

q  =  7.0  X  10“^f0. 07016  +  0.07899X-  0.01358X^  +  0.01757X^ 
s  ‘ 

-  0. 00204X^^4  0.00025  X^  ]  ,  (4.15) 

in  which 

X  =  [3.35  X  10'\,j,  -  1.26  X  10^]  4  400  .  (4.  16) 

The  expressions  for  and  derived  above  are  incorporated  into 
the  w-equation  (2.30)  in  place  of  the  ratio,  H/c^  . 

5.  Method  of  numerical  integration 

It  has  been  a  general  practice  to  test  weather  prediction  models 
by  performing  numerical  integrations  with  observed  initial  data  and 
comparing  the  predictions  with  the  observed  fields.  Although  this  pro¬ 
cedure  is  essential  when  it  is  proposed  to  employ  the  model  in  routine 
forecasting,  it  has  several  disadvantages,  one  of  them  being  the 
necessity  for  a  major  data  collection  and  analysis  program.  Additionally, 
one  cannot  control  the  elements  entering  into  the.  initial  data,  and  con¬ 
sequently  pertinent  factors  relating  to  model  behavior  cannot  be  simply 

isolated. 

For  these  reasons,  the  model  derived  above  is  to  be  applied  to 
a  hypothetical  atmosphere  contained  between  fixed  vertical  walls 
parallel  to  the  zonal  direction  and  extending  to  infinity.  The  inte  - 
gration  region  was  truncated  in  the  zonal  directions  by  assuming  that 


zonal  periodicity  could  be  employed  as  a  boundary  condition  across 
imaginary  meridionally  oriented  walls.  Based  on  the  limitations  of 
computer  core  storage  and  the  linearized  behavior  of  the  planetary 
scale  waves  (see  section  on  linear  analysis),  the  zonal  dimension 
of  the  integration  region  was  fixed  at  10,000  km.  The  zonal  walls 
were  separated  by  4500  km,  which  is  approximately  40  degrees  of 
latitude. 

Within  this  channel  the  initial  field  of  streamfunction  is 
specified  analytically.  From  the  analytic  functions  the  values  of 
streamfunction  were  computed  at  the  grid  points  formed  by  two  sets 
of  line?,  parallel  and  orthogonal  to  the  zonal  walls  respectively. 

The  lines  were  spaced  250  km  apart  and  yielded  an  array  of  19  X  40 
grid  points.  To  this  array  we  added  two  columns  and  two  rows  en¬ 
closing  the  integration  region  for  use  in  specifying  the  boundary  con¬ 
ditions. 

The  differential  equations  of  the  model  fall  into  two  categories, 
two  dimensional  boundary  value  equations  and  tendency  equations. 

The  most  difficult  problem  encountered  is  the  simultaneous  solution  of 
the  omega  and  continuity  equations.  Individually  these  equations  are 
the  familiar  Helmholtz  and  Poisson  equations  respectively.  For  the 
omega  equation  we  used  Dirichlet  conditions  (w  =  0)  on  the  zonal 
boundaries  and  periodicity  on  the  meridional  boundaries.  The 
periodicity  condition  was  duplicated  for  the  continuity  equation,  but 
the  zonal  boundary  condition  was  of  the  Neuman  type  (8x/ay  =  0) 
expressing  the  vanishing  of  the  normal  component  of  the  velocity 

potential  there. 
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The  vorticity  equations  are  of  the  Poisson  type  in  the  stream - 
function  tendency  and  are  solved  for  boundary  conditions  identical  with 
those  for  omega. 

The  tendency  equations  were  of  two  types.  The  streamfunction 
tendency  was  computed  at  each  grid  point  by  the  solution  of  the  vorticity 
equaltion.  The  temporal  extrapolation  is  accomplished  by  first  forward 
and  then  centered  time  steps  over  an  interval  At  =  1  hour.  The  water 
vapor  conveyance  equation  required  the  computation  of  the  Jacobians 
by  means  of  upwind  difference  schemes  to  reduce  the  truncation  error. 
The  tendencies  were  extrapolated  by  forward  differences  only. 

The  finite  difference  equations  used  in  the  numerical  solution 
are  listed  below.  The  notations  are  used  relative  to  this  stencil  of 
grid  points. 


0 


H  +  P3  +  P4  -  4P^ 


"  ^1  -  ^3 


J^(P,  Q)  =  (P^  -  ^3X02  ■  ^4^  ■  ^^2  ■  ■  ^3^ 

J2(P.Q)  =  (P5  -  -  (^8  '  ^7)%  '  ^8’*^4 

J^(P.  Q)  =  (P^  -  P2^^6  ■  ^^4  ■  ^3^^8  ^^2  '  ^3^^7  '  ^^1  '  ^4'^5 

J(P.  Q)  =  j  Q)  +  ^2^^’ 


(5.  1) 
(5.  2) 

(5.3) 

(5.4) 
(5.  5) 
(5.  6) 


=  Qz  -  ^4 

VP  .  VQ  =  [A^P  •  A^Q]  +  [AyP  •  AyQ] 


(5.7) 
(5.  8) 


In  these  equations,  an  unsubscripted  variable  is  to  be  understood 
to  apply  at  the  point  O  in  the  stencil.  The  grid  spacing  will  be  denoted 


by  the  letter  h  (-  250  km). 
The  mean  vorticity  equation: 

W  4 


(5.  9) 


The  shear  vorticity  equation: 


(5.10) 


The  continuity  equation: 


=  -  Vh^w 


(5.11) 


The  omega  equation: 
2 


W 


,  A  h‘ 
2-  f  c 


1 


NEPm.^  EF  ym 


Y_L^j(>)>,v);^))  +(L^j(4.^,^)4j(4>^,^) 


\4h 

E^f.. 


J  J 


->r » ■  A  j  ■  - 


S. 

c 

p 


(5.12) 


In  the  water  vapor  conveyance  equation,  we  make  use  of  upwind 
approximations  for  the  Jacobian  operator  [6].  This  will  be  denoted  by 
the  symbol,  J(4',P)  ■ 

Since  =  +v,  the  northward  component  of  the  wind  is  +v, 

and  =  -u,  the  eastward  component  of  the  wind  is  +u,  we  first 

dy 

determine  the  sign  of  u  and  v  by  using  a  centered  difference  approxi¬ 
mation  for  94j/ax  and  84^/8y  .  We  evaluate  ap/Bx  and  8p/8y  by 
one-sided  differences  in  the  direction  from  which  the  wind  is  blowing. 
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The  water  vapor  conveyance  equation 


8q  1  V  —  _  T  r  ^  _ 

5r  "  ■  TTT  *1)  •  —2  +  (1  -  6)  LE  q  YW 

2h  2h  ^  ^ 


The  auxiliary  equations  are 


r  = 


(5.  13) 


(5.14) 


»T  ^  y 

n 


(5.  15) 


—  =  Q  +  Q 
Cp  s  «c 


1  f  2 

Qs  =  4-81  X  lo-'^r-i^V(^-i|^^)-^F(4^-4.^)l  .(t^-^-^4;  } 


(5.  16) 


(5.17) 


Qg  =  6  17  |w| 


(5.  18) 


The  program  is  written  to  produce  maps  of  the  predicted  fields 
at  6  hour  increments,  as  well  as  a  record  of  the  energy  transformations 
which  are  printed  out  at  hourly  time  steps. 

Because  of  the  inherent  truncation  error  in  numerical  com¬ 
putation,  we  cannot  assure  preservation  of  the  integral  properties  of 
the  model  differential  equations  by  their  numerical  analogues.  We 
have  introduced  the  approximation  to  the  Jacobian  operator  defined  in 
(5.6)  in  an  attempt  to  preserve  the  integral  property  of  the  Jacobian 
operator. 


6.  Linear  analysis  of  the  model  equations 

The  basic  model  equations,  (2.21,  2.24,  2.27,  and  2.29),  may 
be  subjected  to  a  linear  analysis  by  means  of  the  superposition  of  a 
perturbation. 


3-1 


'^rj,  =  gik{x-ct) 

tn  =  iOi 
X'p  ”  i  X ' 


H  =  iH' 


upon  the  basic  state, 


4j  =  -  Uy 
4't  =  -  U,j,  y 
(0  =  0 
Xiji  —  0 

H  =  0 


(6.  1) 


(6.2) 


The  differential  equations  reduce,  upon  substitution  of  (6.1)  and  (6.2), 
to  this  set  of  algebraic  equations  in  the  amplitudes, 


U  + 


U  +  (1 


u. 


y  -  V  -  r  * 


=  0 


(6.3) 

(6.4) 


(C  -  U)+'  +  U.j,$+ =  +  ^  H'  ,  (6.5) 

P 

X'  "  "V  w  =  0  (6.  6) 

It  is  to  be  noted  that  the  term,  •  VXrj.  •  in  equation  (2.  29)  is 
of  higher  order  in  the  perturbation  quantities  and,  consequently,  does 
not  appear  in  (6.  5). 

For  the  present  we  will  examine  the  adiabatic  equations,  taking 

I 
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H'  to  be  zero.  In  order  to  simplify  the  computation  of  the  stability 
characteristics  of  this  system  we  may  first  use  equations  (6.  5)  and 
(6.  6)  to  eliminate  u  and  y'  from  equations  (6,  3)  and  {6.  4).  From 
(6.  5)  and  (6.  6),  we  may  write, 

f  f 

-  if  X'  =  — --2-  [(C  -  U)4>'  +  .  (6.7) 

EFmK 


If  we  set  the  parameter,  which  multiplies  the  bracketed  terms  in  (6.7), 
equal  to  the  dummy  variable,  p'  , 


p*  £  _^  = 


f 

c 

LETmK^. 


(6.  8) 


the  equations  (6.  3)  and  (6.  4)  may  be  rewritten  as 

(c-uA)i.?u,*'  =  o 


and 


Pci  I 

'4^  =  0 


(p'  -  l)U,j,+  +  (1  +  p')(C  -  U)+  (1-  E‘')U^  +-^ 

w  IlC 


The  determinant  of  the  coeficients  of  5  and  i)''  is 


(6.9) 


(6. 10) 


A  =  (1  +  p')X  + 


(1- 


- 


,Pc1 


X  +  (p-  -  l)?u.^^ 


in  which 


X  s 


U  + 


(6.11) 


(6. 12) 


Since  the  system,  (6.9)  and  (6.  10),  is  homogeneous,  a  non¬ 
trivial  solution  exists  only  if  A  vanishes  identically.  Letting  (6.  11) 
equal  zero  and  solving  the  resulting  quadratic  for  the  parameter,  C, 
we  obtain, 


-1 


(6.  H) 


r  ~2 

+  +  (E^  -1)U 

L  K 


2(K^  f  p) 

.  -J 


in  which  the  discriminant,  D,  is  given  by 


D  = 


■pP^  “2 


L  K 


4+  (E  -  1)U, 


-  4 


~2  2 


(6.14) 


If  D  <  0,  the  "phase  velocity",  C,  will  be  complex  and  the 
perturbations  will  be  composed  of  two  linearly  superimposed  parts:  one 
growing  exponentially,  the  other  decaying  exponentially.  The  rate  of 
growth  of  the  amplifying  part  of  the  perturbations  may  be  computed  for 
various  values  cf  the  model  parameters,  as  a  function  of  U,p  and  the 
wave  length,  Zrr/K  . 

Before  discussing  the  graph  showing  the  growth  rate  for  the 
adiabatic  case,  we  may  briefly  consider  the  diabatic  term,  H'.  In  the 
formulation  of  the  diabatic  terms  given  in  section  4,  we  arrived  at  the 
relation 

^  =  IT  |w/  +  Qg  ,  (6, 15) 

P 

in  which  the  first  term  is  the  measure  of  the  release  of  latent  heat 
and  the  second  term  the  rate  of  addition  of  sensible  heat.  The  sensible 
heat  term  is  not  susceptible  to  linear  analysis  since  it  is  a  time 
dependent  non-linear  function  of  the  dependent  variables.  However, 
the  potential  influence  of  the  latent  heat  release  may  be  estimated, 
since  it  is  linearly  related  to  u  when  condensation  is  occurring.  The 
principal  influence  which  we  can  examine  involves  the  effective  increase 
in  the  stability  parameter,  p.  Thus  p  varies  from 
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p  =  0.64  X  lO'^^  cm'^ 

for  the  no -heating  case  to, 

,  ,,  ,--15  -2 

p  =  1.11  X  10  cm 

when  condensation  is  occurring.  By  examination  of  equation  (6,14),  we 
may  determine  the  short-wave  cut-off  to  the  unstable  waves  as  a  function 
of  p.  It  is  clear  that  D  will  be  positive  for  any  value  of  U,p  provided  that 

2 

^<1  .  (6.16) 

Using  the  values  of  p  given  above  we  find  that  for  the  no -heating  case  all 

3 

waves  with  wave  length  less  than  2.49  X  10  km  are  stable.  With  heating, 

the  reduction  in  the  value  of  p  carries  the  potentially  unstable  region  down 

3 

to  waves  with  wave  length  greater  than  1.85  X  10.  km. 

If  we  assume  that  D  is  negative  then  we  may  write, 

C=Cj^±iC^  .  (6.17) 

Using  (6. 13),  it  follows  that 


and 


Pc 

C  =  U  -  4  +  - y- 

^  K  2(K 


/^Pc 


c  ,  .Z2 


+  (£“  -  1)  U„ 


(6.18) 


2(K‘'  +  p) 

If  we  set  equal  to  a  constant,  say  N,  then  we  may  solve 

(6.19)  for  U,j,  as  a  function  of  K.  The  parameter  N  is  the  time  required 
for  wave  number,  K,  to  grow  to  f  -times  its  initial  amplitude  for  the 
value  of  U,j,  as  a  function  of  K  by. 
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in  which 


and 


AC\ 


1/2  1 


2-| 

2.25  -  2.0  E-j 

K^J 


B  = 


pP 


2K 


C  = 


2a  2 

P  P. 


4N 


/'l 


(  6,20) 

(6.21) 

(  6.22) 

(  6.23) 


Equation  (6,20)  was  solved  for  U.j,  ,  with  N  =  g- day,  1  day, 

2  days,  and  10  days,  and  for  the  two  values  of  p  given  above,  was 

-13  -1  -1 

assigned  the  value,  1.45  X  10  .  cm  sec  .  The  results  are  given  in 


figure  1. 


In  addition  to  these  computations  we  also  varied  the  value  of 
— 2  ~2 

the  parameter,  E  ,  We  made  computations  for  E  =  0.333,  0.50,  0.65,  1.00, 
“2 

1.25.  As  E  increased  the  wave  length  of  maximum  instability  shifted 
to  larger  values  and  the  longer  waves  in  general  became  more  unstable 
for  definite  values  of  U,p  . 

The  stability  characteristics  of  the  linearized  model  are  used 
as  a  guide  in  specifying  the  initial  distribution  of  4*  and  in  the 
non-linear  numerical  integrations. 

Two  results  obtained  by  following  the  analysis  procedure  em- 
ployed  by  Wiin-Nielsen  [15]  may  be  noted.  By  solving  the  initial  value 
problem  posed  by  the  linear  ecjuations  for  a  simple  initial  distribution 
of  the  perturbations  I  we  find  that  the  unstable  waves  possess  a  limiting 
phase  relationship  between  the  mean  and  shear  streamfunctions.  If 
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I 


0  represents  the  phase  lag  (positive  when  lags  vjj),  the  limiting 
value  is  given  by 


6  =  arc  tan 


4? 


-  1 


■Pp 

K 


(?  -  1)U„ 


1/2 


(6.24) 


By  evaluating  (6.24)  before  performing  the  numerical  integrations, 
we  are  able  to  set  the  phase  lag  between  4^  and  4^rp  appropriately. 

This  permits  us  to  avoid  computing  changes  which  merely  reflect  a 
readjustment  between  the  two  fields. 

Finally,  use  was  made  of  the  frequency  equation  (6.  13)  to 
compute  the  propagation  velocity  for  the  long,  stable  waves.  It 
was  desired  to  avoid  unrealistic  retrogression  of  the  long  waves  in 
the  computational  scheme  without  introducing  a  mean  divergence  field. 
To  achieve  this  the  zonal  dimension  of  the  integration  region  was  set 
at  10,000  km.  This  value  was  chosen  on  the  basis  of  the  phase  speeds 
for  the  mean  and  shear  streamfunctions  computed  from  (6.  13)  for 
U  =  25  m  sec  ^  and  U,p  =  20  m  sec  ^  (see  Table  1), 


TABLE  1.  The  phase  velocities  of  the  mean  streamfu.iction,  C  , 

and  the  shear  streamfunction,  C-p,  as  a  function  of  wave 
length,  L.  Negative  values  indicate  motion  to  the  west. 


km] 

C  [m  sec  ^  ] 

[m  sec'^  ] 

10 

-  1.66 

+  5.30 

11 

-  9.36 

+  4.47 

12 

-17.79 

+  4.28 

13 

-26.96 

+  4.09 

14 

-36.86 

+  3.99 

15 

-47.49 

+  3.90 

7.  Propoaed  numerical  experiments 

Since  the  energetic  consistency  of  the  model  may  be  tested 
without  the  inclusion  of  the  diabatic  heat  terms,  it  is  proposed  to  carry 
out  comparative  predictions  with  and  without  the  energy  conversion 
term  ( VXij.  *  V'l'rp)  in  the  omega  equation.  For  these  tests,  we  will 
vary  the  initial  stream  function  using  wave  numbers  two  and  three, 
separately  and  in  combination. 

Following  these  adiabatic  computations,  we  will  run 
cases  with  diabatic  heating,  using  both  the  energetically  consistent 
and  inconsistent  versions  of  the  base  model. 

Finally,  we  will  perform  a  set  of  experiments  using  the 
energetically  consistent  version  of  the  model  with  varying  fields  of 
diabatic  heating.  These  last  experiments  are  designed  to  investigate 
the  possibility  that  a  coupling  of  the  diabatic  heating  to  the  adiabatic 
process  of  baroclinic  development  can  augment  the  latter  in  a  sig¬ 
nificant  manner. 
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Fig.  1.  Stability  characteristics  of  the  model  showing  the  wind  shear 
required  for  a  perturbation  of  any  wave  length  to  grow  to 
e  (2.718'««)  times  its  initial  value  in  1/2,  1,  2,  and  10  days. 
Figure  A  is  computed  for  the  parameter,  p,  equal  to  its  value 
in  the  "dry  model",  while  figure  B  is  for  the  p  equal  to  its 
value  in  the  "wet  model"  when  condensation  is  occurring. 
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